#include "sor.h"


void sor(
	double omg,
	double dx,
	double dy,
	int    imax,
	int    jmax,
	double **P,
	double **RS,
	double *res,
	int **Flag,
	double dp
)
{
	int i , j;
	double a , sum = 0;
	int num_of_fluid_cells = 0;

	a = 2.0 * ( (1.0 / pow(dx , 2) ) + ( 1.0 / pow(dy , 2 )));

	/* boundary strip initialization */
	for(j = 1; j <= jmax ; j++) {
		if ( (Flag[0][j] == P_L) && (Flag[imax + 1][j] == P_R) ) {
			P[0][j] = 2*dp - P[1][j] ;
			P[imax + 1][j] = - P[imax][j];
		}
		else {
			P[0][j] = P[1][j];
			P[imax + 1][j] = P[imax][j];
		}
	}


	for(i = 1 ; i <= imax; i++) {
		P[i][0] = P[i][1];
		P[i][jmax + 1] = P[i][jmax];
	}

	/* interior boundary cells initialization */
	for ( i = 1; i <= imax; i++ ) {
		for ( j = 1; j <= jmax; j++ ) {
			switch (Flag[i][j]) {
				case B_N:
					P[i][j] = P[i][j + 1];
					break;
				case B_W:
					P[i][j] = P[i - 1][j];
					break;
				case B_S:
					P[i][j] = P[i][j - 1];
					break;
				case B_E:
					P[i][j] = P[i + 1][j];
					break;
				case B_SW:
					P[i][j] = ( P[i - 1][j] + P[i][j - 1] ) / 2.0;
					break;
				case B_SE:
					P[i][j] = ( P[i + 1][j] + P[i][j - 1] ) / 2.0;
					break;
				case B_NE:
					P[i][j] = ( P[i][j + 1] + P[i + 1][j] ) / 2.0;
					break;
				case B_NW:
					P[i][j] = ( P[i - 1][j] + P[i][j + 1] ) / 2.0;
					break;
				default:
					break;
			}
		}
	}

	for(i = 1 ; i <= imax ; i++) {
		for(j = 1 ; j <= jmax; j++) {
			if ( is_interior(Flag[i][j]) ) {

				P[i][j] = (1.0 - omg)*P[i][j] + (omg/a)*(((P[i+1][j]+P[i-1][j])/pow(dx,2))+((P[i][j+1]+P[i][j-1])/pow(dy,2)) - RS[i][j]);
				num_of_fluid_cells ++;
			}
		}
	}

	for(i = 1 ; i <= imax ; i++) {
		for(j = 1 ; j <= jmax; j++) {
			if ( is_interior(Flag[i][j]) ) {

				sum += pow(((P[i-1][j] - 2*P[i][j] + P[i+1][j])/pow(dx,2))+((P[i-1][j] - 2*P[i][j]+P[i + 1][j])/pow(dx,2)) - RS[i][j],2);
			}
		}
	}

	*res = sqrt(sum/(num_of_fluid_cells));
}
